Chaotic desynchronization of multi-strain diseases 
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Multi-strain diseases are diseases that consist of several strains, or serotypes. The serotypes may 
interact by antibody-dependent enhancement (ADE), in which infection with a single serotype is 
asymptomatic, but infection with a second serotype leads to serious illness accompanied by greater 
infectivity. It has been observed from serotype data of dengue hemorrhagic fever that outbreaks 
of the four serotypes occur asynchronously. Both autonomous and seasonally driven outbreaks 
were studied in a model containing ADE. For sufficiently small ADE, the number of infectives of 
each serotype synchronizes, with outbreaks occurring in phase. When the ADE increases past a 
threshold, the system becomes chaotic, and infectives of each serotype desynchronize. However, 
certain groupings of the primary and secondary infectives remain synchronized even in the chaotic 
regime. 



Recently there has been a wide range of work on 
chaotic synchronization in dynamical systems 0, 
When synchronizing chaotic systems, almost all of the 
work deals with coupled or connected systems and 
analyzing their stability. In biological systems, such as 
population models, synchronization may result from cou- 
pling strengths being enhanced Q, while desynchroniza- 
tion may take place as a result of vaccine control as in 
measles 5J. In this work, we consider the dynamics of 
a single population to shed light on the phase dynamics 
of multi-strain diseases . The dynamics observed ex- 
hibits phase-locked regular behavior, as well as chaotic 
phase desynchronization between strains. Although we 
consider a single population model, we use the term "syn- 
chronization" to describe phase locking between variables 
0. 

Many population models in the past have considered 
single strain diseases, as in childhood diseases. In this 
case, the population may be grouped into the following 
compartments: susceptibles, infectives, and recovered |8j. 
With no seasonal forcing included in the model, the only 
endemic solution to the single strain SIR models is an 
equilibrium point 

However, many diseases have co-circulating strains, or 
seroty pes , such as influenza ^(j; malaria and dengue 
virus |l2(. Such diseases display anti-genic diversity, ex- 
hibiting distinct serotypes when measured. Recent ef- 
forts at modeling multi-strain diseases have explored the 
oscillatory dynamics generated by multiple co-existing 
serotypes with partial cross- immunity How- 
ever, current thinking regarding the interacting serotypes 
of dengue virus is that cross-reactive antibodies act to en- 
hance the infectiousness of a subsequent infection by an- 
other serotype [l5j. This is known as antibody- dependent 
enhancement. 

It has been shown through recent serology measure- 
ments in Thailand that dengue fever, which has four co- 
circulating serotypes, exhibits asynchronous outbreaks. 
That is, each serotype has peaks that occur at different 



times [l^. (See Fig.[7|in the Appendix.) Note that most 
observed infections are secondary [lfjj, due to increased 
symptom severity. 

In this paper, we analyze how the antibody-dependent 
enhancement (ADE) factor controls the onset of oscilla- 
tory outbreaks, as well as how asynchronous secondary 
infections are controlled dynamically. 



I. DESCRIPTION OF MODEL 

To model the spread of multi-strain diseases, we follow 
the approach of Ferguson et. al |l7j , where they restrict 
the model to two serotypes. Our modeling approach dif- 
fers in the general number of serotypes and in that all 
compartments are distinct from one another. The full 
model for n serotypes is described below. Our simula- 
tions will be based on four serotypes, based on measured 
dengue data in Thailand [l6[ . 

The variable definitions are as follows: s, Susceptible 
to all serotypes; Xi, Primary infectious with serotype i; 



Primary recovered from serotype 



1 1 j ■ 



Secondary 



infectious, currently infected with serotype j , but previ- 
ously had i (i ^ j). The model is a system of ODEs 
describing the rates of change of the population fractions 
within each compartment (l8| : 



ds 
~dt 

dxi 
~dt 

drj 
~dt 



0s \Xi + (fi^Xjj -aXi- [i d Xi (2) 



axi - f3n ^2 x j + 4>'^2 x kj ] - Hdn (3) 



2 



dxi 



dt 




HdXij . (4) 



The parameters /x, /i^, (3, and er denote birth, death, con- 
tact, and recovery rates, respectively. We assume that 
individuals who have recovered from two infections are 
immune to further infection since tertiary infections are 
reported very rarely 01 ■ The fixed parameters through- 
out the paper are given by: fj, = fj, d = 0.02, (3 = 200, 
and a = 100, all with units of years -1 (Mortality 
rate, /id, is set equal to the birth rate so that the pop- 
ulation remains constant in time.) Antibody-dependent 
enhancement is governed by the parameter cj>, which has 
not previously been measured for populations. Notice 
that in Eas. 11141 ADE enters in a nonlinear enhancement 
factor when <p > 1. We use a single <j> f° r au strains for 
ease of analysis. Thus any loss of synchrony between the 
strains will result not from asymmetry but from the dy- 
namics itself. Finally, notice that since the value of /x 
is small compared to (3 and a, it can be considered as a 
small parameter. 



II. RESULTS 
A. Bifurcation structure 

Unlike the usual SIR models for single strains, which 
in the absence of forcing have only steady state behav- 
ior, the addition of multiple serotypes can induce regular 
and chaotic outbreaks. In particular, for a critical value 
of <fi, there exists a Hopf bifurcation to periodic oscilla- 
tions. See the bifurcation diagram given in Fig.^for the 
transition from steady state to oscillatory behavior as a 
function of <j>. The usual trivial steady state, which has 
the population consisting of all susceptibles (s = 1) and 
the rest of the components at zero, is unstable. (The 
trivial solution of Eqs. 11141 with n = 4 serotypes has 4 
unstable, 12 strongly attracting, and 5 weakly attract- 
ing directions.) The non-trivial, or endemic, steady state 
may be computed numerically for arbitrary </>. At steady 
state, we notice the following: 1. The primary infectives 
are equal. 2. The recovered variables are equal. 3. All 
secondary infectives are equal. 

Compartmental equality at steady state holds before 
the Hopf bifurcation point as well as after the Hopf point 
(although past the Hopf bifurcation point the steady 
state solution is unstable). We make these assumptions 
about the model at equilibrium, and the resulting local 
dynamics can be reduced to a four-dimensional system: 



dyi 
dt 

dm 

dt 
dy3 
dt 



/j,-nj3y 1 y 2 -n(n-T)l3<py 1 y4, (5) 
2/12/2 + (n-l)/30 yi2/4 -crj/2 (6) 
a y 2 - (n - 1) /3 y 3 y 2 - (n - l) 2 f3 <p y 3 y 4 (7) 




FIG. 1: Bifurcation diagram for the autonomous model with 
(3 = 200, a = 100, ii = 0.02. Shown for each are the 
maxima (black) and minima (gray) of the susceptibles during 
a 100 year run, after transients were removed. 



^ = & 2/32/2 + (n - l)f3(by 3 y i - oy± 



(8) 



Notice that for simplicity we have removed the mortality 
terms in each of the variables, since they are of C(/i) and 
have a negligible effect on the steady states. Moreover, 
removing the mortality terms allows an analytical esti- 
mate of the endemic steady states and stability. Mortal- 
ity does need to be included in the long time asymptotic 
runs, which we do below. The reduced model has the 
following steady state solution: 



P (<p + l)' na' (n-l)/3 (0 + 1)' n(n-l)a 



(9) 



Given the steady state solution as a function of 4> in 
Eq. El to compute the stability we need to evaluate the 
linearization about the steady state. Therefore, we take 
the Jacobian of the vector field of the reduced model 
about the steady state and examine the characteristic 
polynomial for the eigenvalues. Recalling that fi is a 
small parameter, we can expand the solutions to the char- 
acteristic equation in terms of /i. Since the data in [l6| 
displays four serotypes, the number of serotypes is set to 
n = 4. We have a strongly attracting direction given by 
Zi(fi) = —a + O(n), a weakly attracting direction given 
by z 2 (p) = -(3/3 /x l) 2 )/(4cr) + 0{y 2 ), and a pair of 
complex eigenvalues: 



z±(m) 



± z(^) 3/2 (i - fm + eV), (io) 



where //(</>) < 0. The sign of the expression D(<fi) = 
3cj) 2 — 4 <fi — 4 determines the stability of the complex 
pair. Since (j) is assumed to be greater than or equal 
to unity, D(<f>) < if <f> £ [1,2) and positive otherwise. 
Therefore, the steady state undergoes a Hopf bifurcation 
at (j) = 2. The results are close to numerical simulation, 
since mortality terms were dropped in the analysis but 
included in the simulations. 
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FIG. 2: Time series plots of the susceptibles for the au- 
tonomous case, where (3 = 200, a = 100, fj, = 0.02. (a) 
Periodic case for ADE factor <f) — 1.725. (b) Chaotic case for 
4> — 1.73. (c) Longer time series for chaotic case, <f> — 1.73, 
with sampling once a year. 
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FIG. 3: Phase difference analysis of time series in a chaotic 
region with <j> — 1.73 for autonomous system, (a) Time 
intervals between local maxima of secondary infective group 
£2,1. (b) Phase differences between compartments currently 
infected with serotype 1. Primary infective x\ and secondary 
infectives 2:3,1 and £4,1 are compared to 2:2,1 • Maxima occur 
in phase, (c) Comparison between compartments currently 
infected with different serotypes. Secondary infectives 21,2, 
21,3, and 2:1,4 are compared to 2:2,1. Maxima occur out of 
phase during chaotic intervals. 



Since the reduced model does not capture the asyn- 
chronous behavior past the Hopf point, we continue our 
analysis using the full model. As we increase <f> beyond 
the Hopf point, the dynamics exhibits periodic time se- 
ries, as plotted in Fig. Eta). The susceptibles exhibit a 
period of approximately five years when </> — 1.725. The 
actual range of stable periodic solution is quite small, and 
occurs over a A0 of 0.004838. (See the quick transition 
to irregular oscillations after the steady state in the bifur- 
cation diagram given in Fig. Q]) Past the <fi value where 
periodic solutions become unstable, we find chaotic be- 
havior, indicated by a positive maximum Lyapunov ex- 
ponent for most 4> values in that region. (The chaotic at- 
tractors persist over many initial conditions chosen from 
a random distribution.) We note that in this complicated 
region, there are small windows with attracting limit cy- 
cles resulting in a zero maximum Lyapunov exponent. 
Lyapunov exponents were computed by integrating lin- 
ear variational equations along solutions to Eas. 11141 We 
show examples of chaotic oscillations in Fig. [21(b) and (c) 
for short and long time series. Notice that in Fig. |2Jc), 
the time series exhibits oscillatory regions which have a 
slowly growing envelope, interspersed with chaotic inter- 
vals. We will exploit this structure to examine how each 
serotype behaves dynamically. 



B. Phase analysis 

Since the measured data for dengue fever shows that 
the serotypes oscillate out of phase, we investigate the 
phase of primary and secondary infectives with respect 
to a particular secondary infective in the full model of 
Eas. lll4| To measure phase differences with respect to a 
reference infective, let Y(t) denote the reference infective, 
and Z{t) another infective. Let {tk} denote the sequence 
of times for local maxima of Y(t), and {rfc} the times 
for local maxima of Z(t). For r m € (tfe , tfc+i ) , define the 
phase of Z relative to Y in the interval as ^zyijm) — 



In Fig. |31 we compare the relative phases of infective 
groups for <ft = 1.73. For the secondary infective group 
X2,i, we plot the inter-maximum intervals in years in 
Fig. EJa) . Notice that during the non-chaotic times, the 
oscillation intervals grow slowly, until they begin to vary 
in an irregular manner during the chaotic phase. In pan- 
els (b) and (c), a direct comparison between 2:2,1 and the 
other groups is plotted using the phase differencing equa- 
tion ^zY(T m ), normalized between —tt and tt. In panel 
(b), all other infectives who have serotype 1 as the cur- 
rent infection are practically in-phase with group £2.1- In 
contrast, in panel (c), all those having serotype 1 as the 
primary infection, and currently a different serotype as 
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FIG. 4: Dynamics of seasonally driven ADE model, where 
/3o = 200, a — 100, /i = 0.02, and the forcing amplitude 
/3i = 0.05. (a) Periodic time series of susceptible population, 
for the ADE factor <f> =1.5. (b) Quasi-periodic time series 
of susceptibles, for (j> = 1.7244. (c) Projected time series 
sampled at intervals of one year for the susceptible population, 
showing the quasi-periodicity (same parameters as in (b)). 



the secondary infection, lose synchronization when the 
dynamics exhibits chaotic behavior. During the slow 
buildup phase, however, the groups are still synchronized. 
Similar desynchronization during chaotic time series oc- 
curs for the other primary and secondary infectives (not 
pictured here). 



C. Seasonally driven case 

Although the analysis suggests that chaos is respon- 
sible for the observed lag between serotypes, one could 
argue that since there is a seasonal component to the 
disease, adding a periodic forcing term should synchro- 
nize the serotypes, even when they are chaotic. To ad- 
dress this issue, we modified the model to include a con- 
tact rate that modulates with a period of one year; i.e., 
f3(t) = (3q(1 + fli cos(2-7rf)), where j3\ is the forcing am- 
plitude. (/?i = 0.05 was used in this study, but similar 
behavior is observed for other forcing amplitudes.) The 
contact rate prefactor, (3q = 200, is constant as before. 
Analogous to the Hopf bifurcation in the autonomous 
system, bifurcation onto a torus occurs at 4> c — 1.7243. 
For an ADE factor below </j c , we observe periodic be- 
havior, as shown in Fig. 0fa), while for an ADE factor 
just above </j c , we find quasi-periodic behavior, plotted in 
Figs.QIb) and (c). In panel (c), the time series of suscep- 
tibles was sampled at the forcing period and plotted as 
successive iterates to show a cross-section of the torus. In 
both periodic and quasi-periodic cases, the serotypes are 
all in phase, and there is no desynchronization. However, 
for higher ADE, we find that the driven system becomes 
chaotic, and there is desynchronization. 

Figure [S] shows the phase differences between the X2.1 
secondary infective and other infectives, for an ADE fac- 
tor of <f> = 1.74 and forcing amplitude [3\ = 0.05, where 
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FIG. 5: Chaotic phase desynchronization in periodically 
driven system. The ADE factor is cj> = 1.74. (a) Phase differ- 
ences between compartments currently infected with serotype 
1. Primary infective xi and secondary infectives £3,1 and X4,i 
are compared to £2.1- Maxima usually occur in phase, (b) 
Comparison between compartments currently infected with 
different serotypes. Secondary infectives 2:1,2, £1,3, and £1,4 
are compared to 2:2,1. Maxima occur out of phase. (Windows 
of synchrony occur during the oscillatory regions that have a 
slowly growing envelope.) 
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FIG. 6: A histogram plot showing the statistics of the phase 
differences between secondary infections and primary infec- 
tions from Fig.Qj] Black bars: frequency of phase differences 
for compartments currently infected with serotype 1 (data 
from Fig. OJa)), gray bars: frequency of phase differences 
for compartments currently infected with different serotypes 
(data from Fig.EJb)). 



the solution is chaotic. Notice that in the top panel, 
where the phase differences are for other compartments 
currently infected with serotype 1, there is phase syn- 
chrony on average when compared to the case where the 
secondary infections are from a different serotype (second 
panel). Although the phase synchrony is not as good as 
in the autonomous case in Fig. [31 we can get a statistical 
measure showing how on average the phase locking com- 
pares by computing a histogram of both cases. This is 
shown in Fig. where the grey bars correspond to phase 
differences of Fig.JSfb), and the black bars correspond to 
the data from^a). Notice that when comparing primary 
infections of serotype 1 to secondary infections that cur- 
rently have serotype 1, there is a strong phase locking 
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component on average. 

III. CONCLUSIONS 

We have derived and analyzed the dynamics of a model 
for multi-strain diseases with antibody-dependent en- 
hancement. The model for secondary infections, which 
includes ADE as a parameter, adds a new wrinkle to 
models of the SIR type. In previous studies of single 
strain models that do not include environmental forc- 
ing, the endemic equilibrium is the only possible stable 
state. That is, there are no bifurcations which give rise 
to dynamics exhibiting regular or irregular outbreaks. In 
contrast, by modeling the effect of ADE as an increase 
in infectivity of secondary infections, we see both ana- 
lytically and numerically that periodic outbreaks appear 
at a critical ADE value. Moreover, the analysis reveals 
exactly how the period of oscillations depends on the 
ADE parameter near the bifurcation point. The range of 
periods predicted for the parameters used in our compu- 
tations appear to agree well with those observed in the 
data in Fig. in the appendix. 

When the ADE factor increases above a threshold, 
the system's behavior is chaotic, and outbreaks of dif- 
ferent strains occur asynchronously. This observation 
corresponds qualitatively with epidemiological data on 
asynchronous outbreaks of dengue fever (see Appendix). 
Seasonal forcing, thought to be a primary driver for the 
observed oscillations in the different strains, is typically 
believed to disrupt any out-of-phase behavior in the dy- 
namics and force the entire system to lock on the period 
of the forcing. However, in our preliminary study, we 
find that this is not the case. Phase desynchronization 
between serotypes occurs even in the seasonally forced 
case. 

However, there exists a specific relationship between 
the primary and secondary infections. Specifically, we 
have observed that although the different serotypes 
desynchronize when the solutions are chaotic, there 



is surprising structure in the peak outbreaks of the 
serotypes when comparing the appropriate secondary in- 
fectives to the appropriate primary infectives. Although 
there is no vaccine currently available for all serotypes, 
the results here point to potential new methods of anal- 
ysis and monitoring of multi-strain diseases. In the field, 
the majority of the cases reported are secondary infec- 
tions. Therefore, by observing a small percentage of the 
incidence in the secondary infections of one serotype, syn- 
chronization would imply that the data is representative 
of the general behavior of all the groups infected with 
that serotype, including those with only a primary infec- 
tion. Further global analysis techniques based on center 
manifold methods can be used to explain the synchro- 
nization of particular primary and secondary infectives 
when the time series becomes chaotic; this approach is 
the subject of further study. 
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APPENDIX A: EPIDEMIOLOGICAL DATA 

Fig- El reprinted from |l6(, shows the frequency of de- 
tection of each of the four dengue types at one hospital 
in Bangkok, Thailand, over a continuous 27 year period 
of monitoring. Infecting serotypes were defined by iso- 
lation of replication- competent virus and/or detection 
of viral genome in peripheral blood. (It should be noted 
that serological measurements were performed for only a 
fraction of all dengue cases.) Observe that peaks of the 
dengue virus types are asynchronous. 
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FIG. 7: Frequency of detection of each of the four Dengue virus types per month at the Queen Sirikit National Institute for 
Child Health from 1973-1999. Reprinted from [r|. 
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